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Abstract. We comment on an analysis by Contopoulos et al. which demonstrates 
that the governing six-dimensional Einstein equations for the mixmaster space-time 
metric pass the ARS or reduced Painleve test. We note that this is the case irrespective 
of the value, /, of the generating Hamiltonian which is a constant of motion. For J < 
we find numerous closed orbits with two unstable eigenvalues strongly indicating that 
there cannot exist two additional first integrals apart from the Hamiltonian and thus 
that the system, at least for this case, is very likely not integrable. In addition, we 
present numerical evidence that the average Lyapunov exponent nevertheless vanishes. 

The model is thus a very interesting example of a Hamiltonian dynamical system, 
which is likely non-integrable yet passes the reduced Painleve test. 
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1. Introduction. 



In general, it is a difficult task to establish whether or not a Hamiltonian system 
of degrees of freedom is integrable. We note that there are many definitions of 
integrability [|[]. Here, we take as definition [0] the existence of independent ffist 
integrals Ij = Ij{p, q) = const which are in involution, i.e. their mutual Poisson brackets 
vanish {Ii,Ij} = 0, for all i,j. In case that such integrals have not been found, 
one may try to establish the integrability of the dynamical system via the singularity 
analysis called the Painleve test. According to the "Painleve conjecture" a Hamiltonian 
dynamical system is integrable if it has the Painleve property. However, since carrying 
out the Painleve test is often too difficult a task, in practice the test is performed in a 
reduced form called the ARS-algorithm P, ^, ^. In this paper we shall comment on an 
example where the (ARS) conjecture fails. 

In Contopoulos et al. P], and in Cotsakis et al. 0, it is shown that the Einstein 
equations for the mixmaster space-time metric pass the ARS-test (the reduced Painleve 
property) and, by the above proposition, that it is probably integrable. 

Here we shall argue that this is very likely not the case. Our argument pursues 
the following line: First, there is a generating Hamiltonian, the value of which, /, 
characterizes different qualitative behaviours of the model. The case / = corresponds 
to the proper mixmaster universe. Although it contains the chaotic Gauss-map in its 
time-evolution there is no direct contradiction in calling the system integrable, e.g. in the 
sense of being solvable by an inverse scattering transform (see e.g. 0). For / > 0, the 
evolution is in some sense 'trivial', but for / < Of we shall show that the system displays 
chaotic behaviour in the sense that it contains an abundance of unstable periodic orbits, 
even apparently infinite families of these, which in turn renders the existence of a full 
set of ffist integrals very unlikely. In fact, the proliferation of unstable periodic orbits 
indicates the existence of a 'Smale horseshoe' prohibiting integrability. A rigorous proof 
would require finding a transversal homoclinic intersection []T0|, 0, i.e. a transversal 
intersection of the stable and unstable manifolds of a hyperbolic closed orbit of the 
system, a work currently in progress. Second, the Painleve analysis for the mixmaster 
model as presented in Contopoulos et al. does not distinguish the value of this ffist 
integral, i.e. by this analysis one arrives at the contradictory statement that even in the 
chaotic case / < the model is probably integrable:}:. 



f At first sight a universe model with / < might appear quite artificial. It has, however, a physical 
interpretation as a mixmaster universe in which matter with negative energy density is uniformly 
distributed. It gives the possibility of an oscillatory evolution of the three-volume of the universe ^ 
and thus the possibility of periodic (cyclic) mixmaster universes. 

I In a subsequent paper Contopoulos et al. cast doubt as to the conclusion of their original paper 
They introduce a perturbative Painleve test which the mixmaster model fails to pass. Also, Latifi 
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Finally, we present numerical evidence that, despite the existence of unstable 
periodic orbits, the average Lyapunov exponent vanishes in the case / < 0, contradicting 
data based on single trajectories presented in Rugh M and Hobill et al. in . 



2. The equations of motion. 

The mixmaster space-time metric is a famous cosmological toy-model ||T5[ which 
has been studied extensively in the past two decades in various contexts. It is an 
anisotropic generalization of the standard cosmological model of our universe (in case 
our universe is closed), the compact Friedman- Robertson- Walker model. The mixmaster 
metric exhibits very complicated dynamical behaviour near its curvature singularities, 
illustrating an interesting non-linear aspect of the Einstein equations. For a recent 
discussion of the characterization of chaos in general relativity and the mixmaster toy- 
model gravitational collapse, see e.g. S.E.Rugh ||TB|. 
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The governing Einstein equations for the mixmaster metric are given by a set of 
three second order ordinary differential equations 



(1) 



where a, (3 and 7 are the so-called scale factors of the metric, and a dot denotes derivative 
with respect to the (logarithmic) time variable r, see Landau and Lifshitz §118. The 
equations (1) admit a first integral 

I = A{a/3 + + ja) - e^° - e^^ - e^^ + 2(e'("+^) + e2(^+^) + e2(^+")). (2) 

The mixmaster universe, with its interesting oscillatory behaviour of the scale functions 
near its space-time singularity, arises exactly at / = 0. The cosmological model was 
first analyzed by V.A. Belinski and I.M. Khalatnikov and independently by CW. 
Misner [|l^, who showed that the asymptotic evolution could be accurately described 
as a simple combinatorial evolution of the axes governed by the Gauss map (see also 
Barrow and Mayer [^). The appearance of the Gauss map may seem to indicate 
chaotic behaviour of the system but, as the transformation to the relevant coordinates 
is singular and maps non-closed orbits into closed orbits, it is not clear to what extent 
conclusions regarding the inherent chaos of the Gauss map can be carried back to the 
analysis of the differential system governing the mixmaster universe. 

Following the notation of Contopoulos et al. [H, we first perform a change of 



et al. ||l^ note that the model does not pass the full Painleve test. 
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variables : 

X = 2a y = 2f3 z = 2^ , , 

Px = -{y + z) py = -{z + x) pz = -{x + y) 

and find that the system can also be described as governed by the Hamiltonian 

H = \{pl+pl+pl-2p.,Py-2pyp,-2p,p,)W''W^+e^-2e^^y-2ey+'-2e'+^ = -/.(4) 

For the purpose of the Painleve analysis a further change of variables is made 

X = e^ Y = ey Z = e^ (5) 
and we finally end up with this version of differential equations: 

2X = X{p^-Py-p^) 2Y = Y{py-p^-p^) 2Z = Z{p^-p^-py) 
p^ = 2X{Y + Z-X) py = 2Y{Z + X -Y) p, = 2Z{X + Y - Z), 

in which the first integral takes the form 



(6) 



I = -lipl+Pl + Pl) + \iPx+Py+P.r - 2(^^ + Y' + Z') + {X + Y + Zf. (7) 

The sign of the first integral characterizes three different phases of the model. 
This can be seen from the scale invariance of the differential equations (j^) under the 
transformation 



(8) 



In popular terms we may say that / = characterizes a phase transition for the 
system, as the quahtative description of the time evolution depends strongly on the sign 
of J IQ. In particular we note here that since the time variable scales with |/|~^/^, the 
Lyapunov exponent for the system will have to vanish for J — > 0, revealing an apparent 
non-chaotic behaviour at / = 0, i.e. for the mixmaster universe. Due to this scale 
invariance we are free to choose an arbitrary value of / when we look for periodic orbits 
of the mixmaster system for / < 0, and we will choose the value / = — 1. 
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3. Painleve analysis. 

The Painleve analysis, when applied to (H), showed that the system passes the (ARS) 
Painleve test for the case / = and therefore, according to the cited conjecture, that 
the system is probably integrable. The analysis does not, however, depend on the value 
of /. Here we focus on one of the singular solutions (cf. p. 5798 in Contopoulos et al. 
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Inserting and identifying coefficients in the Laurent series expansion show that the 
singular expansions 

X = ^+7iS+... Y = X2S + ... Z = X3S + ... 

Px = + b + CS + . . . Py = P2 + 2iX2S + . . . Pz= P3 + 2iX3S + . . . 

satisfy the set of equations (||) to zero'th order in s = t — tq, provided b = p2 + Ps and 
c = 2i[x2 + X3 — 271]. Here, t is the complex time variable and Tq, 71, X2, ^3,^2 and ps 
are the six free parameters in the generic Painleve expansion. Computing / to zero'th 
order (all other terms must vanish due to the time invariance of /) we get : 

= + 2i(x2 + X3 - 371) , (10) 

which is a constant of motion but in general does not vanish. Thus we see, as is also the 
case for the other classes of singular solutions in 0, that the set of equations has the 
(ARS) Painleve property independent of the value of the ffist integral /. 



4. Closed orbits and non-integrability. 

One way to characterize non-integrability is by means of non-zero Lyapunov exponents. 
In general, each ffist integral of a Hamiltonian system implies the vanishing of two 
Lyapunov exponents due to preservation of symplectic two- forms ||2T|. In practice this 



characterization is not very useful. First, it is not easy to establish numerically due 
to the infinite time limit involved, and second, because the reverse need not be true, 
i.e. the vanishing of 2n Lyapunov exponents does not imply the existence of n first 
integrals. However, in general the vanishing of Lyapunov exponents occurs for almost 
all trajectories of infinite length. In particular, we shall consider closed orbits, i.e. 
trajectories which return to a given point, x, in the phase space after some period 
T. Periodic orbits have the great advantage (see also [^]) that they are completely 
determined by their finite time behaviour and therefore provide a numerically easy and 
reliable way of stating non-integrability of the system. 

We write the differential equations (^-(^ in the form : 

= 1^ ° ^ VH{x) = JVH{x) , X = {p,q) e R'^ X (11) 

(here d = ?>) and we denote by (j)H{x) the flow obtained by integrating this equation. A 
real-analytic function / on i?^'' is said to be a ffist integral if it commutes with H, i.e. 
if their poisson bracket vanishes : 

{I,H} = {VI,JVH) = (12) 

(■, ■) being the usual inner product in R^"^. In particular 0, this condition implies that 
the flow 0f (x), generated by JV/, commutes with 0^. Thus, we have the two relations : 
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/ o 0*^ = / te R 

0J o 0*^ = 0*^ o s,teR 



(13) 



Assume now that we have found a periodic orbit x = (pfjix) of our Hamiltonian flow 
with return time T > 0. We may use this periodicity condition in our two relations 
above, taking derivatives in x and in s (at s = 0), respectively, to obtain (tr meaning 
transpose) : 

VlixfDf^ix) = Vlixf , 
JV/(x) = D(P'^{x)JVI{x) . 

We can replace / with H or any other first integral (provided they all commute) and 
conclude that the (transposed) set of VH, V/, etc. generates a set of left invariant 
eigenvectors of D(pjj{x) while by multiplying by J we obtain a set of right invariant 
eigenvectors of the same matrix. 

The two sets are evidently mutually orthogonal (as the first integrals commute). 
Thus we may orthonormalize each set and take their union to be a basis of R^'^. A 
straightforward calculation shows that in this basis D(j)fj{x) takes the form of an upper 
(or lower) triangular matrix with 1 in the diagonal. Thus, we conclude that if the system 
has d independent first integrals at a periodic orbit then the stability matrix has only 
one eigenvalue, 1, with degeneracy 2N. 

As Poincare already noted , the presence of non-unit eigenvalues does not exclude 



integrability but could germinate from a degeneracy of first integrals. The simplest 
example of this is the stiff pendulum which, when placed at rest in its top position, is 
in an unstable fixed point. The separatrices through this point separate the phase 
space into oscillatory modes and two types of rotational modes. The eigenvalues of the 
stability matrix for the fixed point for any non-zero time are not on the unit-circle. The 
arguments above do not apply to this case as we here have VH = (a special case of 
degenerating first integrals). More complicated situations arise in integrable systems in 
higher dimensions when, for closed orbits, the gradients of the first integrals become 
linearly dependent. 

As mentioned above, a rigorous proof would require the presence of a transversal 



homoclinic intersection [1^, |TT|. However, we insist that an abundance of unstable 
periodic orbits strongly indicates non-integrability of the dynamical system, i.e. that 
it cannot possess d commuting analytic first integrals which are independent almost 
everywhere. 
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5. Unstable periodic orbits for / < 0. 

Finding periodic orbits in a low-dimensional dynamical system such as (1) means solving 
the equation 

= X (15) 

and is a straightforward, if tedious, task. One starts with random initial conditions 
and lets the system evolve until a suitable close-return occurs. What is a suitable close- 
return is as much determined by elapsed time as distance between initial and final values 
of the variables, due to the high instability of the dynamics. One then uses the initial 
conditions of the close return trajectory as a starting guess in an extended Newton 
search for a solution to ([T5|). This takes the following steps: 

1. Choose a Poincare section, i.e. a {2N — l)-dimensional surface that the trajectory 
will cross. This is rather arbitrary, but a plane through the initial point perpendicular 
to the initial velocity works well. 

2. Integrate the equations of motion including the Jacobian matrix M (here v = JVH) 

M = ^^M, M,.„ = 1, (16) 

until the trajectory crosses the Poincare section at xj, close to the initial point, Xj. 

3. Assume that the trajectory is close to a periodic orbit and solve the linearized version 
of (USD: 



X = (f)^{x) ^ 0^(xi) + M{x - Xi) =Xf + M{x - Xi), (17) 

i.e. solve 

{1- M)x = Xf - Mxi. (18) 

This presents the problem that since M has 1 as (at least) a double eigenvalue, 1 — M 
is not invertible. To remedy this, we augment the equation (^) with 2 more equations 
and 2 more variables The two extra equations serve the purpose of keeping the 



solution X on the Poincare section and preserving the first integral / during the Newton 
iteration, while the two variables add corrections to x along the right/left eigenvectors 
of the eigenvalue 1, i.e. v{x) and V/, in order to insure the existence of a solution to 

4. Repeat until satisfactory convergence is obtained. 

We plot two examples of periodic orbits in figures |I] and ||. In tables |l] and ^ we 
give the initial conditions, periods and eigenvalues of a few of the periodic orbits found 
through this method. As can be seen, the eigenvalue 1 has multiplicity 2 in all our 
examples. We have found approximately 200 unstable periodic orbits, of which around 
50 were found in a random search through phase space, whereas the remaining number 
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was found in a systematic search for members of infinite families of orbits as described 
below. This numerical result makes it highly unlikely that the system is integrable for 
J < 0. 

The existence of unstable periodic orbits does not directly lead to chaos in the sense 
of a positive Lyapunov exponent in the mixmaster system for / < 0. In fact, we have 
been able to find series of periodic orbits that indicate that the Lyapunov exponent will 
converge to zero for infinite time. Consider a family of periodic orbits where one of the 
scale factors, a, makes one oscillation while the other two, P and 7 make n oscillations. 
In figure § such an orbit with n = 10 is shown. We find that for such a family the 
largest eigenvalue of M, A„ grows linearly with n, whereas the period, T„ grows as y/n. 
This means that for the instability exponent A„, i.e. the local Lyapunov exponent of 
the individual periodic orbits, we have: 

_ logA^ logn 



which will tend to zero for n going to infinity. This is shown in table |^ and figure ^. The 
limit of such a family of periodic orbits is therefore marginally stable. In the computation 
of the Lyapunov exponent, marginally stable regions seem to dominate making the 
exponent converge to zero. This conclusion is supported by a direct computation of 
the finite time approximation of the Lyapunov exponent. Starting with 1000 random 
initial conditions, we have computed the average finite time Lyapunov exponent, which 
appears to fall off as a power law with increasing r. This result is shown in figure |^. 

We therefore conclude that the mixmaster system (1) for /< is neither integrable 
nor chaotic in the sense of having a positive maximal Lyapunov exponent. 

6. Concluding remarks 

We have noted that, according to the recent analysis by Contopoulos et al. 0, the 
mixmaster cosmological model passes the Painleve test not only for the case / = 0, i.e. 
the vacuum Einstein equations for which the behaviour of the axes is described by the 
BKL-combinatorial model, but for arbitrary values of the first integral /. At the same 
time we have shown that when / < 0, the model is probably non-integrable through the 
existence of an abundance of unstable periodic orbits. 

Thus, the model provides a first example of a Hamiltonian system of non-linear 
differential equations that passes the ARS-test for the Painleve property but which 
is not integrable. One may still speculate whether the mixmaster model has the full 
Painleve property (the answer appears to be negative according to the recent results of 
Contopoulos et al. and those of Latifi et al. |jl3[), or whether the constraint 1 = 
may be incorporated in a dynamical context or in the Painleve analysis in some way, 
avoiding the non-integrable periodic solutions described in this article and rendering the 
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system integrable. 
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Table 1. Initial conditions for some periodic orbits. Entries i, v and vi correspond 
to figures 0, H+js] and IJ. The number of digits given corresponds to the accuracy of 
the numerical methods employed as estimated by using 2 different integration routines 
and 2 different tolerance levels in the routines. 







Initial values 




Period 


# 


a, a 




7, 7 


T 


i 


-0.7465064136 


-0.04110593028 


-0.31459136109 


5.6040350442 




0.0 


0.57358291296 


-0.6330334642 




ii 


0.28088191655 


-1.4043414748 


-0.21644617646 


7.9472870699 




0.0 


-0.19132448256 


0.084543111864 




iii 


-0.85822654286 


0.296817610585 


-0.374109406546 


8.54414096178 




0.0 


0.0 


0.0 




iv 


0.399801735 


-0.141272123 


-0.26896997 


13.4135432128 




0.0 


1.044336138 


-0.47263976 




V 


0.13651016517 


-0.83335865226 


-0.5987346358 


14.6263998708 




0.0 


0.424450201779 


-0.32431285874 




vi 


0.67215626848 


-0.04851348437 


-0.04851348437 


16.15030271 




0.0 


0.2302667762 


-0.230266776 




vii 


-1.0663696240 


-2.7978105880 


0.07780069395 


16.435444052 




0.0 


0.0668561349 


0.3463999832 





Table 2. Eigenvalues for the periodic orbits of Table 0. The lesser accuracy of the 
unit eigenvalues is caused by the fact that they appear as double eigenvalues. 



Eigenvalues 


i 


-29.654128172 






1.00000 


-0.948098014 + 


i0.3179782317 




-0.033722117682 




1.00000 


-0.948098014 - 


iO. 3179782317 


ii 


-103.545671685 






1.00000 


-0.69755560723 + 


iO. 71653065169 




-9.657574128 10' 


-3 


1.00000 


-0.69755560723 - 


iO.71653065169 


iii 


136.451552742 






1.00000 


-0.48706495012 + 


iO. 87336575062 




7.328608432 10" 


-3 


1.00000 


-0.48706495012 - 


i0.87336575062 


iv 


2540.8960825 






1.00000 


1.044452134 






3.9356194 


10- 


-4 


1.00000 


0.957439758 




V 


-5843.077705 






1.0000 


0.5307158319 + 


iO.8475498249 




-1.7114268 


10- 


-4 


1.0000 


0.5307158319 - 


iO.8475498249 


vi 


-643.8934046 






1.00000 


0.0021357407 + 


iO. 9999977193 




-1.55305209 


10- 


-3 


1.00000 


0.0021357407 - 


i0.9999977193 


vii 


1926.43285564 






1.0000 


2.170605581 






5.1909413 


10- 


-4 


1.0000 


0.4607009253 





Table 3. Largest stability eigenvalue, period and stability exponent for the {l,n,n) 
family of periodic orbits. T„ and A„ vs. n are plotted in figure ^. Beyond n w 40 it 
gets increasingly difficult to follow the trajectories numerically because the dynamics 
is determined by the decreasing difference between the two oscillating scale factors P 



and 


7- 








n 


T 


A 


A 


^ iog(|A|)/r 


1 


5.6040350441 


-29.654128172 


-0, 


,11430453548 


2 


7.9472870699 


-103.5456716852 


-0, 


.04481126745499 


3 


9.42569727483 


-186.757808548 


-0, 


.0280031805718 


5 


11.7015648632 


-331.02446237 


-0, 


.0175279864076 


10 


16.150302715 


-643.89340465 


-0, 


.0100444159622 


15 


19.649946461 


-967.1141240 


-0, 


.0071080716705 


20 


22.618966753 


-1292.2369183 


-0, 


.0055439756742 


25 


25.242560127 


-1617.915587 


-0, 


.00456692177 


30 


27.61864459 


-1943.8400780 


-0, 


.00389559861571 


40 


31.84403793 


-2596.030034 


-0, 


.003028369681 


46 


34.12932414 


-2987.456574 


-0, 


.002678592127 
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Figure Captions 



Figure 1. The simplest periodic orbit of the mixmaster system corresponding to entry 
i in table The 3 scale factors a, f3 and 7 each make one identical oscillation. The 
trajectory takes the form of a slightly distorted circle. 

Figure 2. A more complex periodic orbit of the mixmaster system corresponding to 
entry v in table |l|. The largest eigenvalue A of the Jacobian for this trajectory is w 
5800 for a period T of only 14.62. 

Figure 3. The same trajectory as above but shown projected on the (a,/?)- 
plane. 

Figure 4. An example of a periodic orbit in the family with respectivelt 1, n and n 
oscillations of the three scale factors with n = 10. The initial values correspond to 
entry vi in table |l|. For increasing n the difference between the two oscillating scale 
factors becomes smaller and smaller. 

Figure 5. The power law dependence of the largest instability eigenvalue A of the 
Jacobian, M, and the period, T, vs. n in the {l,n,n) family of periodic orbits. We 
have found orbits up to n — 46. For n > 6 we estimate the slope to be .495 for T and 
1.006 for A. Some of the data points are listed in table ^j. 

Figure 6. The finite time Lyapunov exponent Ar computed as an average of 1000 
random trajectories. For numerical reasons we use the largest eigenvalue of M^^M. 
r = 70 is approximately the time limit, where we lose information about the lowest 
eigenvalue due to numerical underflow. We believe this explains the bend of the curve 
near this value and therefore that the power law fall-off holds even for r > 70. The 
power law fall-off indicates that the inflnite time Lyapunov exponent is 0. 
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